%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Dispersive light atom interaction with Schrödinger & Bloch.
%  n  = (n11 n22 n12 n21) .
%
%  ---	 |3>
%    \
%    --- |2>
%    /
%  ---	 |1>
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function Dispersive
    figure(1);					clf;
    V=[0:.1:2*pi];

    subplot(221),plot3(sin(V),cos(V),0*V,'b','MarkerSize',10,'LineWidth',1.2)
    set(gca,'xLim',[-1 1],'yLim',[-1 1],'zLim',[-1 1],'FontSize',18,'LineWidth',1.2)
    hold on
    subplot(221),plot3(0*V,sin(V),cos(V),'b','MarkerSize',10,'LineWidth',1.2)
    set(gca,'xLim',[-1 1],'yLim',[-1 1],'zLim',[-1 1],'FontSize',18,'LineWidth',1.2)
    subplot(221),plot3(cos(V),0*V,sin(V),'b','MarkerSize',10,'LineWidth',1.2)
    set(gca,'xLim',[-1 1],'yLim',[-1 1],'zLim',[-1 1],'FontSize',18,'LineWidth',1.2)
    
    G12=  .01;
    G23=    1;
    G13=    1;
    O12=    1;
    O23= 1000;
    d12=    0;
    Dt1=  .04;      Dt2= .001;
    T1 =   39;      T2 =   39;
    D23=O23^2/4/pi*Dt2*T2;

    Schrod(G12,G23,G13,O12,O23,d12,D23,Dt1,Dt2,T1,T2);
    Bloch(G12,G23,G13,O12,O23,d12,D23,Dt1,Dt2,T1,T2);
end

function Schrod(G12,G23,G13,O12,O23,d12,D23,Dt1,Dt2,T1,T2)
    psi0=[1 0 0]';
    col='rg';
    for nn=0:1
        psi=psi0;
        for tt=1:(2*T1+T2)
                o12=O12;        dt = Dt1;
                o23=  0;        d23=1e-9;
            if tt<(T1+T2)
                o12=  0;        dt = Dt2;
                o23=O23*nn;     d23= D23;
            end
            if tt<T1
                o12=O12;        dt = Dt1;
                o23=  0;        d23=1e-9;
            end

            psi = Hamilton_K(psi,dt,o12,o23,d12,d23);
            psi = Hamilton_K_Pert(psi,dt,o12,o23,d12,d23);

            n(1)=-2*imag(psi(1)*psi(2)');
            n(2)=-2*real(psi(1)*psi(2)');
            n(3)=real(abs(psi(2)^2)-abs(psi(1)^2));

            subplot(221),plot3(n(1),n(2),n(3),['.',col(nn+1)],'MarkerSize',10,'LineWidth',1.8)
            set(gca,'xLim',[-1 1],'yLim',[-1 1],'zLim',[-1 1],'FontSize',18,'LineWidth',1.2)
    %         title(['|\rho| =  ',num2str(norm(n)),'|\rho| =  ',num2str(norm(rho))])
            xlabel('\fontsize{18}Im {\it\rho_{12}}')
            ylabel('\fontsize{18}Re {\it\rho_{12}}')
            zlabel('\fontsize{18}2{\it\rho_{22}}-1')
            text(-1.5,-1,1.5,'\fontsize{16}(a)')
            view([30 40])
            pause(.01)
        end
    end
end

function Bloch(G12,G23,G13,O12,O23,d12,D23,Dt1,Dt2,T1,T2);
    V=[0:.1:2*pi];

    subplot(222),plot3(sin(V),cos(V),0*V,'b','MarkerSize',10,'LineWidth',1.2)
    set(gca,'xLim',[-1 1],'yLim',[-1 1],'zLim',[-1 1],'FontSize',18,'LineWidth',1.2)
    hold on
    subplot(222),plot3(0*V,sin(V),cos(V),'b','MarkerSize',10,'LineWidth',1.2)
    set(gca,'xLim',[-1 1],'yLim',[-1 1],'zLim',[-1 1],'FontSize',18,'LineWidth',1.2)
    subplot(222),plot3(cos(V),0*V,sin(V),'b','MarkerSize',10,'LineWidth',1.2)
    set(gca,'xLim',[-1 1],'yLim',[-1 1],'zLim',[-1 1],'FontSize',18,'LineWidth',1.2)
    
    rho0=[1 0 0 0 0 0 0 0 0]';
    col='rg';
    for nn=0:1
        rho=rho0;
        for tt=1:(2*T1+T2)
            o12=O12;            dt = Dt1;
            o23=  0;            d23=   0;
            if tt<(T1+T2)
                o12=  0;        dt = Dt2;
                o23=O23*nn;     d23= D23;
            end
            if tt<T1
                o12=O12;        dt = Dt1;
                o23=  0;        d23=   0;
            end

            rho = Matrix_K(rho,dt,G12,G23,G13,o12,o23,d12,d23);

            n(1) =2*imag(rho(4));
            n(2) =2*real(rho(4));
            n(3) =real(rho(2)-rho(1));

            subplot(222),plot3(n(1),n(2),n(3),['.',col(nn+1)],'MarkerSize',10,'LineWidth',1.8)
            set(gca,'xLim',[-1 1],'yLim',[-1 1],'zLim',[-1 1],'FontSize',18,'LineWidth',1.2)
    %         title(['|\rho| =  ',num2str(norm(n)),'|\rho| =  ',num2str(norm(rho))])
            xlabel('\fontsize{18}Im {\it\rho_{12}}')
            ylabel('\fontsize{18}Re {\it\rho_{12}}')
            zlabel('\fontsize{18}2{\it\rho_{22}}-1')
            text(-1.5,-1,1.5,'\fontsize{16}(b)')
            view([30 40])
            pause(.01)
        end
    end
end

function psi = Hamilton_K_Pert(psi0,dt,o12,o23,d12,d23)
% Schrödinger equation with perturbation Hamiltonian
    H=[
        0   	o12/2               0 
        o12/2   -o23^2/4/d23-d12    0
        0		0                   d23+o23^2/4/d23
    ];

    psi = expm(1i*H*dt)*psi0;
end

function psi = Hamilton_K(psi0,dt,o12,o23,d12,d23)
% Schrödinger equation with exact Hamiltonian
    H=[
        0         o12/2       o12/2
        o12/2     0           o23/2
        0		  o23/2       d23
    ];

    psi = expm(1i*H*dt)*psi0;
end

function n = Matrix_K(n0,dt,G12,G23,G13,o12,o23,d12,d23)
% Bloch equation
    I=          1i/2;
    g12=         G12;                                                       % decay of coherences
    g13=     G13+G23;
    g23= G12+G13+G23;
    d13=     d12+d23;
    L12=g12/2+1i*d12;
    L13=g13/2+1i*d13;
    L23=g23/2+1i*d23;
    
    M=[
        0		G12    G13    I*o12   -I*o12	0		0		0       0
        0	   -G12    G23   -I*o12	   I*o12	0		0		I*o23	-I*o23
        0		0	  -G13-G23  0		0		0		0		-I*o23	I*o23
        I*o12  -I*o12	0   	-L12	0		I*o23	0		0		0
        -I*o12  I*o12	0       0		-L12'	0		-I*o23	0		0
        0		0		0		I*o23	0		-L13	0		-I*o12	0
        0		0		0		0		-I*o23	0		-L13'	0		I*o12
        0       I*o23  -I*o23	0		0		-I*o12	0		-L23	0
        0      -I*o23	I*o23	0		0		0		I*o12	0		-L23'
    ];

    n = expm(M*dt)*n0;
end
